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AN ENTROPY PRESERVING FINITE-ELEMENT/FINITE- VOLUME PRESSURE 
CORRECTION SCHEME FOR THE DRIFT-FLUX MODEL 

L. Gastaldo 1 , R. Herbin 2 and J.-C. Latche 3 

Abstract. We present in this paper a pressure correction scheme for the drift-flux model combining 
finite element and finite volume discretizations, which is shown to enjoy essential stability features 
of the continuous problem: the scheme is conservative, the unknowns are kept within their physical 
bounds and, in the homogeneous case (i.e. when the drift velocity vanishes), the discrete entropy of 
the system decreases; in addition, when using for the drift velocity a closure law which takes the form 
of a Darcy-like relation, the drift term becomes dissipative. Finally, the present algorithm preserves 
a constant pressure and a constant velocity through moving interfaces between phases. To ensure the 
stability as well as to obtain this latter property, a key ingredient is to couple the mass balance and 
the transport equation for the dispersed phase in an original pressure correction step. The existence 
of a solution to each step of the algorithm is proven; in particular, the existence of a solution to the 
pressure correction step is derived as a consequence of a more general existence result for discrete 
problems associated to the drift-flux model. Numerical tests show a near-first-order convergence rate 
for the scheme, both in time and space, and confirm its stability. 
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1. Introduction 

Dispersed two-phase flows and, in particular, bubbly flows are widely encountered in industrial applications 
as, for instance, nuclear safety studies, which are the context of the present work. Within the rather large 
panel of models dealing with such flows, the simplest is the so-called drift-flux model, which consists in balance 
equations for an equivalent continuum representing both the gaseous and the liquid phase. For isothermal 
flows, this approach leads to a system of three balance equations, namely the overall mass, the gas mass and 
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the momentum balance, which reads: 



dp 

i H 

dpy 
dt 
dpu 
~~dt~ 



- V • {pu) = 

+ V • (py u) = -V • {py (1 - y) u r ) + V • (DVy) 
+ V • (pu <g> u) + X7p - V • t(u) = f v 



(1) 



where t stands for the time, p, u and p are the (average) density, velocity and pressure in the flow and y stands 
for the gas mass fraction. The diffusion coefficient D represents in most applications small scale perturbations 
of the flow due to the presence of the dispersed phase, sometimes called "diphasic turbulence" and u r is the 
relative velocity between the liquid and the gaseous phase (the so-called drift velocity); for both these quantities, 
a phenomenologic relation must be supplied. The forcing term /„ may represent, for instance, the gravity forces. 
The tensor r is the viscous part of the stress tensor, given by the following expression: 



t(u) = p (Vw + V*u) - - p, (V • u) I 



For a constant viscosity, this relation yields: 



V ■ T 



Am 



VV • u 



(2) 



(3) 



and, in this case, this term is dissipative (i.e. for any regular velocity field u vanishing on the boundary, the 
integral of V • t(u) ■ u over the computational domain is non-negative). 

This system must be complemented by an equation of state, which takes the general form: 



Q P ' a (p, a g ) = (1 - OLg)p t + a g Qg{jp) 



(4) 



where a g stands for the void fraction and g g (p) expresses the gas density as a function of the pressure; in the 
ideal gas approximation and for an isothermal flow, g g (-) is simply a linear function: 



Qg(p) 



(5) 



where a is a constant characteristic of the gas, equal to the sound velocity in an isothermal (monophasic) flow. 
The density of the liquid phase pe is assumed to be constant. Introducing the mass gas fraction y in (JH) by 
using the relation a g g g = py leads to the following equation of state: 



(p, y) 



Qg{p)Pl 



>y + (i-y) Qg{p) 



(6) 



The problem is supposed to be posed over fi, an open bounded connected subset of M. d , d < 3, and over a 
finite time interval (0,T). It must be supplemented by suitable boundary conditions, and initial conditions for 
p, u and y. 

To design a numerical scheme for the solution of the system ([lj, one is faced with several difficulties. First, 
since the fluid density pi is supposed not to depend on the pressure, almost incompressible zones, i.e. zones 
where the void fraction is low, may coexist in the flow with compressible zones, i.e. zones where the void fraction 
remains significant. This feature makes the problem particularly difficult to solve from a numerical point of 
view, because the employed numerical scheme will have to cope with a wide range of Mach numbers, starting 
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from zero to, let us say, for low to moderate speed flows, a fraction of unity. Second, the gas mass fraction y can 
be expected, both for physical and mathematical reasons, to remain in the [0, 1] interval, and it appears strongly 
desirable that the numerical scheme reproduces this behaviour at the discrete level. Finally, it appears from 
numerical experiments that, in order to avoid numerical instabilities, the algorithm should preserve a constant 
pressure through moving interfaces between phases (i.e. contact discontinuities of the underlying hyperbolic 
system). To obtain a scheme stable in the low Mach number limit, the solution that we adopt here is to use 
an algorithm inspired from the incompressible flow numerics, namely from the class of finite element pressure 
correction methods, and which degenerates to a classical projection scheme when the fluid density is constant. 
The last two requirements are met thanks to an original pressure correction step in which the mass balance 
equation is solved simultaneously with a part of the gas mass balance. For technical reasons, the solution of 
this latter equation is itself split in two steps, the first step thus being incorporated to the pressure correction 
step and the second one being performed independently. 

This work takes benefit of ideas developped in a wide literature, so we are only able to quote here some 
references, the choice of which will unfortunately probably appear somewhat arbitrary. For a description of 
projection schemes for incompressible flow, see e.g. [19, 24] and references herein. An extension to barotropic 
Navier-Stokes equations close to the scheme developped here can be found in [14], together with references to 
(a large number of) related works (see e.g. [21] for the seminal work and [30] for a comprehensive introduction). 
Extensions of pressure correction algorithms for multi-phase flows are more scarce, and seem to be restricted to 
iterative algorithms, often similar in spirit to the usual SIMPLE algorithm for incompressible flows [22,25,29]. 
The gas mass balance equation, i. e. the second equation of ((T]) , is a convection-diffusion equation which differs 
from the usual mass balance for chemical species in compressible multi-component flows studied by Larrouturou 
[23] by the addition of a non-linear term of the form V • p<p(y) u r , where <p(-) is a regular function such that 
<p(0) = ip(l) = (in the present case, ip(y) — y (1 — y)). In [17], we propose a finite-volume scheme for the 
numerical approximation of this type of equation, and we prove the existence and uniqueness of the solution, 
together with the fact that it remains within physical bounds, i.e. within the interval [0, 1]. Here, the proof of 
the same results combines arguments from both [23] and [17]. 

Several theoretical issues concerning the proposed scheme are studied in this paper. First, the existence 
of a solution to the pressure correction step, which consists in an algebraic non-linear system, is obtained by 
a topological degree argument. Second, we address the stability of the scheme. At the continuous level, the 
existence of an entropy for the system when the drift velocity vanishes (i.e. the homogeneous model) is well- 
known. In addition, it is shown in [20], by a Chapman-Enskog expansion technique, that the two-fluid model 
can be reduced to the drift-flux model when a strong coupling of both phases is assumed, with a Darcy-like 
closure relation for the drift velocity, i. e. an expression of the form: 



1 /, \ Q g (p)-Pe 
- (l-a g )a g 

A p 



u r = - (l-a 3 )a 3 ^ Vp (7) 



where A is a positive phenomenological coefficient. The same relation can also be obtained by neglecting in the 
two-fluid model the difference of acceleration between both phases [28]. With such an expression for u r , the 
drift term becomes a second order term, and it is shown in [20] that it is consistent with the entropy of the 
homogeneous model (i.e. that it generates a non- negative dissipation of the entropy). These results are proven 
here at the discrete level: up to a minor modification of the proposed scheme, which seems useless in practice, 
the entropy is conserved when u r is equal to zero, and when the closure relation ([7]) applies and with a specific 
discretization, the drift term generates a dissipation. 

This paper is built as follows. The fractional step algorithm for the solution of the whole problem is first 
presented in section [2] together with some of its properties: the existence of a solution to each step of the 
algorithm, the fact that the unknowns are kept within their physical bounds and that the algorithm is able 
to preserve a constant pressure and a constant velocity through moving interfaces between phases. The proof 
of the existence of the solution to the pressure correction step is obtained as a consequence of a more general 
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existence theory for some discrete problems associated to the drift-flux model, which is exposed in the appendix. 
Next two sections are devoted to the stability analysis of the scheme; after establishing estimates for the work of 
the pressure forces (section [3]), we first address the case u r = (section [4. ip . then the case where u r is given by 
the Darcy-like closure relation ([7]) (section [4. 2|) . Finally, numerical tests are reported in section [51 they include 
a problem exhibiting an analytical solution which allows to assess convergence properties of the discretization, 
a sloshing transient in a cavity, and the evolution of a bubble column. 

For the sake of simplicity, we suppose for the presentation of the scheme and its analysis (sections [H |3] and 
that the velocity is prescribed to zero on the whole boundary dtt of the computational domain, and that the gas 
mass flux through dQ, so both u r and the component of normal to the boundary, also vanishes. Moreover, 
the analysis of the scheme assumes that pure liquid zones does not exist in the flow, which, with the proposed 
algorithm, is a consequence that such zones are not present at the initial time (i.e., at t = 0, y G (0, 1]); getting 
rid of this latter limitation at the theoretical level seems indeed to be a difficult task. However, the numerical 
tests presented in section [5] are not restricted to theses situations. In particular, y — in the liquid column in 
the sloshing problem, up to spurious phases mixing by the numerical diffusion near the free surface; it is also 
the case at the initial time in the bubble column simulation. 

In the presentation of the scheme, the drift velocity is supposed to be known, i.e. to be given by a closure 
relation independent of the unknowns of the problem, and this still holds in numerical experiments. The case 
where u r is given by ([7]) is thus only treated from a theoretical point of view in section 14.21 



2. The numerical algorithm 



We present in this section the numerical scheme considered in this paper. We begin by describing the 
proposed scheme in the time semi-discrete setting, then we introduce the spatial discretization spaces and we 
detail the discrete approximation and the properties for each step of the algorithm at hand. 

2.1. Time semi-discrete formulation 

Let us consider a partition = to < ti < . . . < tjv = T of the time interval (0,T), which is supposed uniform 
for the sake of simplicity. Let St be the constant time step St — t n+ \ — t n for n = 0, 1, . . . , N. In a time 
semi-discrete setting, the algorithm proposed in this paper is the following three steps scheme: 
1 - solve for u n+1 



p n-l u n 



St 



+ v • ( P n u n ® u n+l ) + v P n - v • r(u n+1 ) = 



(8) 



2 - solve for p n+1 , u n+1 , p n+1 and z n+1 



St 

g p ' z (p n+1 , z n+1 )-p r ' 



V(p n+ -P n ) = o 



St 



V • (g(p 



■" + ' z n+1 )u n+1 ) = 



z n+l _ p n y n 



St 



+ V-(z n+1 u n+1 )=0 



g p ' z (p n+ \z n+1 ) 



(9) 



3 - solve for y n+1 

„n+l n+l n+1 

" + V • {p n+1 y n+1 {\ - y n+1 ) < +1 ) = V • [DVy n+1 ) (10) 
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The first step consists in a classical semi-implicit solution of the momentum balance equation to obtain a 
predicted velocity. 

Step 2 is an original nonlinear pressure correction step, which couples the mass balance equation (second 
equation) with the transport terms of the gas mass balance equation (third equation). A new unknown is 
introduced in this step instead of the gas mass fraction, the partial gas density z given by z = py. Thus, the 
equation of state must be reformulated to express the mixture density as a function of the partial gas density 
and of the pressure, which, from equation ([6]), yields: 



When the liquid and the gas densities are very different, this law presents much less steep variations than 
the relation linking the density and the mass fraction y, specially in the neighbourhood of y — 0; this change 
of variable thus makes the resolution of this step much easier, and the overall algorithm more robust. In 
counterpart, it leads to split the gas mass balance equation: transport terms are dealt with in the present 
step, and the gas mass fraction is corrected in a next step (step 3) to take into account the drift terms. The 
pressure correction step would degenerate in the usual projection step as used in incompressible flows solvers 
if the density was constant (i.e. z — 0). Taking (at the algebraic level, see section |2T4|) the divergence of the 
first relation of ([9]) and using the second one to eliminate the unknown velocity u n+1 yields a non- linear elliptic 
problem for the pressure. Solving at the same time this elliptic problem and the third equation by a Newton's 
algorithm, we obtain the pressure and the gas mass fraction. Once the pressure is computed, the first relation 
yields the updated velocity and the fourth one gives the end-of-step density. 

Finally, in the third step, the remaining terms of the gas mass balance are considered, and the end-of-step 
gas mass fraction is computed. 

The motivations of this time discretization are the following ones: to keep the mass fraction y in the phys- 
ical range [0, 1], to allow the transport of phases interfaces without generating spurious pressure and velocity 
variations and to ensure the stability of (i.e. the conservation of the entropy by) the scheme. To show how this 
time splitting algorithm achieves these goals is the aim of the remaining of this paper. 

2.2. Spatial discretization 

Let M. be a decomposition of the domain either into convex quadrilaterals (d — 2) or hexahedrons (d = 3) 
or in simplices. By £ and £(K) we denote the set of all (d — l)-edges a of the mesh and of the element K £ M. 
respectively. The set of edges included in the boundary of f2 is denoted by £ ex t and the set of internal ones 
(i.e. £ \ £ cx t) is denoted by £- ln t. The decomposition M. is supposed to be regular in the usual sense of the 
finite element literature (e.g. [6]), and, in particular, M. satisfies the following properties: f2 = Uks-M if 
K, L £ A4, then K<1 L = or Kf] L is a common edge of K and L, which is denoted by K \L. For each internal 
edge of the mesh a = K\L, ukl stands for the normal vector of a, oriented from K to L. By \K\ and \a\ we 
denote the measure, respectively, of K and of the edge a. 

For stability reasons, the spatial discretization must preferably be based on pairs of velocity and pressure 
approximation spaces satisfying the so-called inf-sup or Babuska-Brezzi condition (e.g. [4]). Among these 
elements, nonconforming approximations with degrees of freedom for the velocity located at the center of the 
faces seem to be well suited to a coupling with a finite volume treatment of the other equations, as is proposed 
hereafter for the gas mass balance; this is the choice made here. The spatial discretization thus relies either 
on the so-called "rotated bilinear element" /Pq introduced by Rannacher and Turek [26] for quadrilateral or 
hexahedric meshes, or on the Crouzeix-Raviart element (see [7] for the seminal paper and, for instance, [10, p. 
83-85] for a synthetic presentation) for simplicial meshes. The reference element K for the rotated bilinear 
element is the unit d-cube (with edges parallel to the coordinate axes); the discrete functional space on K is 




(11) 
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Qi(K) d , where Q\(K) is defined as follows: 

Qi(K) = span {f , (^)i=i,...,d, (a;? -»f +1 )i=i,... > d_i} 

The reference element for the Crouzeix-Raviart is the unit (i-simplex and the discrete functional space is the 
space Pi of affine polynomials. For both velocity elements used here, the degrees of freedom are determined by 
the following set of nodal functionals: 

{F^i, ae£(K), i = l,...,d}, F^{v) = l^- 1 [ Vidr? (12) 

J a 

The mapping from the reference element to the actual one is, for the Rannacher-Turek element, the standard 
Qi mapping and, for the Crouzeix-Raviart element, the standard affine mapping. Finally, in both cases, the 
continuity of the average value of discrete velocities (i.e., for a discrete velocity field v, F^^v), 1 < i < d) 
across each face of the mesh is required, thus the discrete space Wh is defined as follows: 

W h = {v h e L 2 (n) : v h \ K G Q x {K) d , MK G M; 

F ai i(vh) continuous across each edge a G £ m t, for 1 < i < d; 
F a ,i(v h ) = 0, Vaefcxt, l<i<d} 

For both Rannacher-Turek and Crouzeix-Raviart discretizations, the pressure is approximated by the space Lh 
of piecewise constant functions: 

Lh = {qh 6 L 2 (VL) : q h \ K = constant, MK G M} 

Since only the continuity of the integral over each edge of the mesh is imposed, the velocities are discontinuous 
through each edge; the discretization is thus nonconforming in H 1 (Vl) d . These pairs of approximation spaces 
for the velocity and the pressure are inf-sup stable, in the usual sense for "piecewise H 1 " discrete velocities, i.e. 
there exists q > independent of the mesh such that: 

p V • v dx 

Vp G L h , sup > q ||p - m(p)|| L 2 (n) 

vew h IM|l,6 

where m(p) is the mean value of p over f2, the symbol / stands for N / and stands for the broken 

Jn,h kgm Jk 

\Wib = E / l^| 2 d*= / \Vvfdx 



Sobolev H 1 semi-norm: 



From the definition p^|) . each velocity degree of freedom can be univoquely associated to an element edge. 
Hence, the velocity degrees of freedom may be indexed by the number of the component and the associated 
edge, and the set of velocity degrees of freedom reads: 

{v„ t i, a G £- mt , 1 <i <d} 

We define v a = J2i=i v <r,i e ^ where eW is the i th vector of the canonical basis of R d . We denote by tp^ the 
vector shape function associated to v a ^, which, by the definition of the considered finite elements, reads: 
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where tp a is a scalar function. 

Each degree of freedom for the pressure is associated to a mesh K, and the set of pressure degrees of freedom 
is denoted by {pk, K £ M}. As the pressure, the density p, the gas mass fraction y and the gas partial density 
z are approximated by piecewise constant functions over each element, and the associated sets of degrees of 
freedom are denoted by {px, K € M}, {vk, K € M} and {zk-, K G M} respectively. 

2.3. Spatial discretization of the momentum balance equation 

The main difficulty in the discretization of the momentum balance equation is to build a discrete convection 
operator which enjoy the discrete analogue of the kinetic energy relation, that is: 



n l 



dp u 



Of 



+ V • (pu(g>u) 



u dx 



dt 



1 



■p\u\ 



provided that 



dp 

at 



+ V- {pu) = o 



To this purpose, we follow an idea developped in [2], and already exploited for the same problem as here in [17]. 
The idea is to derive a finite- volume- like discretization of the convection operator, in order to apply the following 
result [14]. 

Theorem 2.1 (Stability of a finite volume advection operator). Let {p* K )xeM an d {pk)k^m be two families 
of positive real numbers satisfying the following set of equation: 



MK g M, 



\K\ 
St 



(PK ~ Pk) + 



E 

a=K\L 







(13) 



where F a _K is a quantity associated to the edge a and to the control volume K ; we suppose that, for any internal 
edge a = K\L, F„ t K — —F<t,l- Let {s* k )k£M an d {sx)KeM be two families of real numbers. The following 
stability property holds: 



E 

KeM 



ZK 



-JfipKSK ~ P* K S* K ) + ]T F a 5 
a=K\L 



SK + S L 



K 



>1 y \K\ 
~ 2 ^ St 



2 * * 2 

PKS K ~ p K S K 



(14) 



K 



a 6 £{K) Y 

a e £(K) 

Figure 1. Diamond-cells for the Crouzeix-Raviart and Rannacher-Turek element. 

To this purpose, we first define a control volume for each degree of freedom of the velocity, that is, in view 
of the discretization used here, around each barycenter of an internal edge. Let a — K\L and Dk.it be the 
conic volume having a for basis and the mass center of K as additional vertex (see figure [1} . The volume 
F> a = F>K.a U Dl,<t is referred to as the "diamond cell" associated to cr and Dk,g is the half-diamond cell 
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associated to a and K. For the Crouzeix-Raviart element and, for the Rannacher-Turek element, when the 
mesh is a rectangle (in two dimensions) or a cuboid (in three dimensions), the integral of the shape function 
associated to the edge a over the element K is the measure of the half-diamond cell Dx a . Thus, the application 
of the mass lumping to the terms of the form pu leads, in the equations associated to the velocity on the edge a, 
to a discrete expression of the form p a u ai where p a results from an average of the values taken by the density 
in the two elements adjacent to a, weighted by the measure of the half-diamonds: 

Vae£ int , \D a \ Pa = \D K , a \ PK + \D L , a \ p L (15) 

where \D a \ is the measure of the diamond cell D a , \Dk.<t\ and |-Dl.<t| are the measure of the half-diamond cells 
associated respectively to a and K and to a and L. This lumped time derivative term naturally combines with 
a discretization of the advective term of the form: 

V<7 € £int, advection term ~ F £ . a u e 

where £(D a ) is the set of the edges of D ai u E is a centered approximation of nonce £{D a ) and F £ _ <7 is a mass 
flux through e. To proceed, we must now derive for this latter quantity an approximation which satisfies the 
compatibility condition p3[) of theorem l2.ll (in fact, the discrete mass balance over the diamond cells). Suppose 
that we are able to build, for any control volume K, a field pu K (x) such that V • pu K (x) remains constant inside 
the element K and that we take for the mass flux F e ^ a through each diamond cell edge e included in K: 

F e ,a = J PUk( x ) ' n e,<? drf(%) 

where n e . a is the normal vector to e outward D a . As the divergence of pu K is constant over K, it may be 
checked that, if the flux of pu K through each edge of K is the same as the mass flux used in a discrete mass 
balance over K, let say |c| (p u) a ■ ?V, this mass balance "carry over" the half-diamond cells Dk,o, which, by 
summation over the two half-diamond cells, yields a compatibility condition of the desired form [17]. Such a field 
pu K (x) is derived for the Crouzeix-Raviart element by direct interpolation (i.e. using the standard expansion 
of the Crouzeix-Raviart elements) of the quantities ((p u)a-)ae£(K) : 

<reS(K) 

For the Rannacher-Turek element, when the mesh is a rectangle or a cuboid, it is obtained by the following 
interpolation formula: 

P^k( x ) = X! oi a (x-n a ) [(pu) a -n a ] n a 

cr££(K) 

where the a a (-) are affine interpolation functions which are determined in such a way that the desired conditions 
hold, i.e. that the flux of pu through each edge a of K is |<r| (p u) a ■ n a . Extension to more general grids is 
underway. Finally, since, in the proposed fractional step algorithm, the mass balance equation is considered 
only when the solution of the momentum balance is achieved, to obtain the desired compatibility condition (|13D , 
we use the mass balance at the previous time step: the approximations of the density in the time derivative 
term are shifted of one time step and the quantities ((p u)(r)creS(K) used to compute the mass fluxes F™ a are 
chosen to be the mass fluxes obtained in the discrete mass balance at the previous time step. Since standard 
finite elements techniques are used to discretize the term Vp™ — V • r(u n+1 ), this yields the following discrete 
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momentum balance equation: 

V<7 € £int, for 1 < i < d, 

l~ Pff l ( n u n+1 - u n ■) + - F n (u n+1 + u 11 ^ 1 ) 

ee£(D CT ), 

where, the bilinear form ad(-, •) represents the viscous term and, Vi> £E W^, Vw £ VF^, is defined as 

follows: 



a d (v,w) 



Vv : Vw + — V • v V • w 

■ ) 



da; if Q holds (case of constant viscosity), 

/ t(v) : \7w dx with r given by © otherwise. 



Note that, for Crouzeix-Raviart elements, a combined finite volume/finite element method similar to the 
technique employed here has already been analysed for a transient non-linear convection-diffusion equation by 
Feistauer and co-workers [1,9,13]. 

As a consequence of the stability of the convection operator, we have the following regularity result. 

Lemma 2.2 (Properties of the numerical scheme - velocity prediction). Let us assume that the viscous term is 
dissipative (i.e. \fv € Wu, o,d(v, v) > 0, which holds for the form of a<i(-, •) used in case of a constant viscosity); 
then the first step of the scheme, namely the velocity prediction step, has a unique solution. 

Remark 2.3 (First time step). To ensure the compatibility condition (fl"3"|) at first time step, a prediction step 
must be used to initialize the density: 

^^+V-(A- 1 ) = (17) 
where p~ x and it -1 are suitable approximations for the initial density and the velocity, respectively. 

2.4. Spatial discretization of the pressure correction step 

The discretization of the first equation of the pressure correction step is consistent with the momentum bal- 
ance one, i.e. we use a mass lumping technique for the unsteady term and a standard finite element formulation 
for the gradient of the pressure increment: 

Va e & mt , for i < i < d, p » (<+' - <Y) + I (p n+1 - p n ) V • <p® dx = 

ot Jn,h 

As the pressure is piecewise constant, the transposed of the discrete gradient operator takes the form of the finite 
volume standard discretization of the divergence based on the finite element mesh, thus the previous relation 
can been rewritten as follows: 

Va££ int , a = K\L, |^ pi « +1 - < +1 ) + \a\ [(p n K +1 - pi) - (p n L +1 - pi)] n KL = (18) 
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Similarly, as the density is piecewise constant, the approximation of the time derivative of the density in the 
mass balance will also look as a finite volume term. This point suggests a finite volume discretization of this 
latter equation, which reads: 

VK e M, J|l (q*'Wc +1 , z n K +1 ) - P \) + J2 Kg = (19) 

a=K\L 

To ensure the positivity of the density, we use an upwinding technique for the convection term, then the mass 
flux from K across a — K\L, F™~^}, is expressed as follows: 

Kg = M (p n+1 u n+ % • n a = (vt, K ) n+1 e p ' z (p n K + \z n K +1 ) - (y- K ) n+1 q^{pT\zI +1 ) 

where (v+ K )™ +1 and (v~ K ) n+i stands respectively for max(v™~^, 0) and — min(v™"^ L , 0) with v™^ = |er| ■ 
n K L- 

Consistently with the mass balance equation, we use for the discretization of the third relation of ([9]), i.e. 
the transport of the gas partial density z, a finite volume method with an upwind technique for the convection 
term V • (zu). This yields the following discrete equation: 

\/K £ M, ^ (z n+1 - P \ yl) + J2 KkT +1 " K. K ) n+1 = (20) 

a=K\L 

In the following lemma, we state some properties of this pressure correction step which are obtained as a 
particular case of the existence theory presented in section El 

Lemma 2.4 (Properties of the numerical scheme - pressure correction step). Let the density of the liquid phase 
be constant and the gas phase obeys the ideal gas law. Then, under the assumption that, WK £ M, p\ > and 
y% £ (0, 1], the system (|18[) - (|20p has a solution, and any solution of this step is such that: 

n+l 

VKeM, p n K +1 > 0, p™ +1 > 0, zl +1 > and £ (0, 1] 

Pk 



Let us now turn to the practical solution of this pressure correction step. Keeping the same notation for the 
unknown functions and the vectors gathering their degrees of freedom, the algebraic formulation of this step 
reads: 

jM p n (u n+1 - u' l+1 ) + B* {p n+1 - p n ) = 

^R(^(p"+\ z n+1 ) - p n ) - BQJ +I u n+1 = (21) 

^R (z n+1 - p n y n ) - BqJ +1 u n+1 = 

In the first relation, M p ™ stands for the diagonal mass matrix weighted by the density at t n (at edge center) p™, 
so the diagonal entry of M p ™ associated to the internal edge a and the component i reads (M^n)^ = \D a \ p™. 
The matrix B* of M. NxM , where N — d card (£mt) and M = card (.M), is associated to the gradient operator; 
consequently, the matrix B is associated to the opposite of the divergence operator. In the second and in the 
third relation, Q u ^ +1 (with w = p or w = z) is a diagonal matrice, the entry of which corresponding to an edge 
a £ £j n t, a = K\L, is obtained by just taking w at t n+1 in the element located upstream of a with respect to 
u n+1 , i.e. either or u>2 +1 . The matrix R is diagonal and, for any K £ A4, its entry R^ is the measure of 

the element K. 
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The elliptic problem for the pressure is obtained by multiplying the first relation of (j!H]) by B Q"^+i (M^n+i ) 1 
and using the second one. This equation reads: 

L p n+1 + ^ R e p ' z (p n+1 . y n+1 ) = l p" + -^ R p n + ^ B Qp- +1 ( 22 ) 

where, as seen in [14], L = B Q u £ +1 (M p n) _1 B* can be equivalently evaluated in the "finite volume way" by the 
following relation, valid for each element K: 

a=K\L P ° 

where p up .<r stands for the upwind density associated to the edge a. One recognize in this relation a usual 
finite volume diffusion operator, with a particular diffusion coefficient which, for instance, can be evaluated for 
rectangular parallelepipedic control volumes as d Pup, a /Pa- The factor d should be suppressed to be consistent 
with what would be obtained by a finite volume discretization of this elliptic equation, if this latter was derived 
in the time semi-discrete setting: this fact is linked with the well-known non-consistency of the Rannacher-Turek 
or Crouzeix-Raviart discretization of the Darcy problem. 

Then, equation (|22|) is solved at the same time as the third equation of l[2"Tj) by a Newton's algorithm. Once 
p^fj is known, the first relation of l|2ip gives the updated value of the velocity: 



U k+1 = U 



^(M^B^^-p-) (23) 



As, to preserve the positivity of the density, we want to use in the mass balance the value of the density 
upwinded with respect to u n+1 , equations (|22[) and (|23[) are not decoupled, by contrast with what happens in 
usual projection methods. They are thus solved in sequence, performing the upwinding with respect to in 
and then updating the velocity by (|23p . up to convergence. 



2.5. Spatial discretization of the correction step for y 

To be consistent with the discretization of the first part of the gas mass balance, the correction step for y is 
discretized by the finite volume method, and the resulting discrete problem reads: 

\k\ 53 (g:+ k Y g(yl + \yl +1 ) - (G? K Y 9(y2 + \y n K +1 ) 

a=K\L ( 24 ) 

d 

=K\L 



In this relation, for all edge a = K\L, d a is the Euclidean distance between two points xk and xj, of the 
adjacent meshes K and L, supposed to be such that the segment [xk, Xl] is perpendicular to K\L. These 
points may be defined as follows: if the control volume K is a rectangle or a cuboid, xk is the barycenter of K; 
if the control volume if is a simplex, xk is the circumcenter of the vertices of K. Note that, in this latter case, 
the condition xk € K implies some geometrical constraints for K . Of course, in the cases where the diffusion 
coefficient D = 0, these limitations are useless. 

The quantities {G n +^)+ and {G n +^)~ are defined as {G n +^)+ = max(G^+ 1 , 0) and (G' CT l + 1 )- = - min{G%£, 0) 
respectively, with G™" 1 ^ given by: 

(~<n+l _ n+1 / n+1 
u a,K — Ptr,up / u r n <J 
J a 
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where p"„p stands for p™+ P = p^ +1 if u™ +1 ■ n G > and p"™ = p2 +1 otherwise. Note that this upwind 
choice with respect to u n+1 has no theoretical justification: in fact, the developments of this paper hold with 
any discretization for this density, and we use here the same discretization as in the mass balance simply to 
make the informatic implementation easier. The function •) corresponds to an approximation of tp(y) = 
max[y (1 — y), 0] by a monotone numerical flux function. Let us recall the definition of this latter notion [12]: 

Definition 2.5 (Monotone numerical flux function). Let the function g(-, •) 6 C(R 2 , R) satisfy the following 
assumptions: 

(1) g(di , 02) is non-decreasing with respect to a\ and non-increasing with respect to 02 , for any real numbers 
d\ and a,2, 

(2) <?(•, •) is Lipschitz continuous with respect to both variables over R, 

(3) g(di, a±) = tp(a,i), for any a\ G R. 

Then <?(•,•) is said to be a monotone numerical flux function for ip(-). 

Several choices are possible for the numerical flux function g(-, •) and we refer to [12] for some examples and 
references. We adopt here the following simple flux-splitting formula: 

g(a 1 ,a 2 ) = gi{a{) + g 2 (a 2 ) 

where gi[a\) — a\ if a\ € [0,1] and gi(ai) — otherwise, and (72(^2) = — i a 2) 2 if 0,2 € [0,1] and 32(^2) = 
otherwise. Note that this choice does not exactly match the definition, as neither gi(-) nor 32(0 are continuous 
at ai = 1. However, this is unimportant, as one can prove, even in this case, that the solution y remains in the 
interval (0, 1], as stated in the following lemma which is a weaker version of the result proven in [17, section 2]. 

Lemma 2.6 (Existence and uniqueness for a discrete solution). Let us suppose that, WK e M, > 

and z^ 1 1 p 7 ^ 1 € (0,1]. Then, there exists a unique solution to the considered discrete problem (|24[) , and this 
solution verifies y^ 1 € (0, 1], VK £ M. 

2.6. Some properties of the scheme 

The following theorem gathers some properties of the scheme, which are essentially straightforward conse- 
quences of lemmas 12.21 12.41 and 12.61 

Theorem 2.7 (Properties of the scheme). Let the density of the liquid phase be constant and the gas phase obeys 
the ideal gas law. We suppose that the viscous term is dissipative (i.e. Vt> € Wh, ad(v,v) >0). Ln addition, we 
assume that the initial density is positive and the initial gas mass fraction belongs to the interval (0,1]. Then 
there exists a solution (u n )i< n <N> (p n )i<n<N , {p n )i<n<N> (z n )i<n<N an d {y n )i<n<N to the scheme which 
enjoys the following properties, for all n < N: 

• the unknowns lie in their physical range: 

VK e M, p n K >0, z%>0, P £>0, 1/5-6(0,1] 

• the total mass, the gas mass and, if f v = 0, the integral of the momentum are conserved: 

E 1*1 Pk ^ E M Pk 

KeM KEM 

E \K\ z n K = E 1*1 PkVk - E 1*1 PkVk 

KeM KEM KEM 

E \D a \ P n - l K= E \ D '\p^< 
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We now turn to another feature of the scheme, which, from numerical experiments, seems to be crucial for the 
robustness of the algorithm. Let us suppose until the end of this section that the drift velocity u r , the diffusive 
coefficient D and the forcing term f v are set to zero. In addition, we make abstraction of the boundary condition, 
i.e. we momentarily reason as if the problem was posed in M™. Then the continuous problem enjoys the following 
property: if the initial velocity and the initial pressure are constant, let say u = u and p = Po respectively, then 
they remain constant throughout the transient, while p or z are transported by this (constant) velocity; this 
solution corresponds to the transport of the contact discontinuity of the underlying hyperbolic system, the wave 
structure of which is quite similar to the Euler equations one [20]. The objective of the subsequent development 
is to prove that the numerical scheme considered in this paper presents the same behaviour: if, at the initial 
time, u° K = uq and p\ — po for all K E Ai, then p^ 1 = po and u^ 1 — uo, for all K S M. and n < N. 



Let us assume that, at time t — t n , the velocity u n and the pressure p n take the constant value uo and po 
respectively. We are now going to check that there exists a solution u n+1 , p n+1 , z n+1 and y n+1 to the scheme 
such that u n+1 — uq and p n+1 — pq. The discrete momentum balance equation reads, with a zero forcing term: 



Ver G & ln t, for 1 < i < d, 



ee£(D a 
e=D a \D, 



Replacing u n and p n by uq and po respectively and taking , u™ +1 = uq for all a G £i n t, this system becomes: 



Vct g £ int , a = K\L, 



"o 



\D* 

St 



see{D a ) 



which is verified thanks to the equivalence between mass balances over primal and dual meshes, as explained in 
section |2~31 We now turn to the pressure correction step, which we recall: 



\D 



^p n a « +1 -K +1 ) + H [(pk +1 -pk)-(pI +1 -pI)] n KL = o. 



St 

~St ^' </ ' ,x ' ' h 



Q^{p n K + \z n K +1 )-p n K ] 



\K 



(z n+1 



E [«kY +1 4 +1 - ( 



v - ) n + 1 z n + x 



Vcr G £ int , o = K\L 

VK eM 
VK eM 



cre£(K) 



Taking = u n for all a E £ m t and p 7 ^ 1 = po for all K E M, the left hand side of the first equation of this 
system vanishes. Next, following [15], we remark that, at fixed pressure, the equation of state giving the density 
p as a function of z becomes an affine function: 

p=gP>z(p ,z) = z(l-^) +p e 
V PO J 
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Introducing this relation in the mass balance equation, we obtain: 



\K\ 
St 



> a 



+ E «k 

<y££(K) 



PO 
vn+1 



y n+l 
-K 



Pi a 
Pa 



+ Pi 



n+1 



v n+l 



Pa 



which can be recast as: 



1 - 



Pi a 
Pa 



\K\ 
St 



(Zk +1 -Zk)+ E [«K) n+1 Zk +1 ~ K, K ) n+1 



<re£(K) 



+pz E <i 1 = 

<re£(K) 



which, as the last term vanishes for u n+1 = uq, is exactly the same equation as the gas mass balance. Thus, 
u n+i _ p n+1 = p 0j z n+1 given by this latter equation and y n+1 satisfying the correction step (which, 
for u r = and D = becomes p n + 1 y n + 1 = z n+1 ) is a solution to the scheme. Consequently, provided that 
the solution is unique, the algorithm does preserve constant pressure and velocity through moving interfaces 
between phases, and transport this interface with this constant velocity. 

Remark 2.8 (More general boundary conditions). The same property holds with a bounded computational 
domain when prescribing on the boundary either u = uq or a Neumann condition compatible with u = uo and 
p = po; this fact has been confirmed by numerical experiments, although we leave its proof beyond the scope of 
this presentation, to avoid the technicalities of the description of these latter discrete boundary conditions. 

Remark 2.9 (On the choice of coupling the mass balance and the gas mass balance equations). As in [17], 
one may be tempted, specially for computing efficiency reasons, to use a fully fractional step algorithm, i.e. to 
solve all the equations sequentially. The central argument of the preceding development is that, with a fixed 
pressure, the quantity py is affine with respect to p, and this fact originates from the particular form of the 
equation of state. Thus, for this argument to hold, it is mandatory for the density in the product py to be given 
by the equation of state p = Q p ' v {y,p*), where only the pressure may be taken at the previous time step or at 
the previous stage of the algorithm (indeed, when checking as below that the interface is transported with a 
fixed pressure, p will be considered constant, in particular with respect to time). Hence, the transport terms in 
the gas mass balance should read, in the time semi-discrete setting: 

1 (Q™(y n+1 ,p n ) y n+1 - P n y n ) + V • Q™(y n+1 ,p n ) y n+1 u n + • • • = 

But, in this case, the compatibility condition which yields a maximum principle for the advection operator, 
which here would read: 

j t (e p ' v (y n+ \p n ) - P n ) + v • e ™{y n+ \p n )u n = o 

does not hold. So it seems that an algorithm keeping y within its physical bounds and transporting the interface 
at constant pressure and velocity necessarily couples the mixture and the gas mass balance. 



3. The stability induced by the pressure forces work 

The aim of this section is to prove that the discretization at hand satisfies a stability bound which can be 
seen as the discrete analogue of the following equation : 

- / p(x) V • u{x) dx=^- [ f(x) dx (25) 
Jn dt J n 
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where f(x) stands for the volumetric free energy of the mixture. The role played by this estimate in the theory 
which is developped here is twofold. First, it provides an a priori bound for a class of discrete problems including 
the pressure correction step, which is the corner stone to prove the existence of a solution; this development is 
presented in appendix. Second, it is crucial to derive stability results for the scheme. 

Throughout this section, we suppose that both the drift velocity and the gass fraction diffusion vanishes, so 
the overall and gas mass balance equations simply read: 

(26) 
(27) 

Of course, stability results for the complete problem will in fine depend on the fact that the neglected terms in 
(|2"Tj) are dissipative with respect to the free energy; this point will be treated further. 

This section is organized as follows. First, we prove this estimate in a general setting, i.e. without specifying 
the equation of state for the fluid. Then we explain as this theory applies to the case specifically adressed here, 
namely a constant density fluid and a gaseous phase obeying the ideal gas law. 



dp 
dt 
dz 
dt 



V ■ (pu) 

V ■ (zu) 



3.1. Abstract estimates 



The formal computation which allows to derive estimate ([23)) in the continuous setting is the following. The 
first assumption is that, through the (system of) equation(s) of state, the specific free energy can be expressed as 
a function of the mixture density and the gas partial density, which we write / = f(p,z). Then multiplying the 
mass balance equation by the derivative of / with respect to p, the gas mass balance equation by the derivative 
of /(■, •) with respect to z and finally summing these relations, we obtain: 



df 
dp 



df 

dz 



= 



which yields: 



d df df 

Q- t f{p{x, t), z(x, t)) + -± V • (pu) + -± V • (zu) = 



Developping the divergence terms, we get: 



g^f(p(x,t),z(x,t)) + u- 



df „ df „ 

dp dz 



V • u 



d A+ z d X 
dp dz 



= 



(28) 



The second term of this relation is equal to u ■ V f(p(x,t), z(x,t)). Adding and substracting / V • u, we thus 
have: 



d_ 

di 



f(p(x, t), z(x, t)) + V • (f(p(x, t),z(x, t)) u)+V-u 



df , df , 

P ^~ + z T~ ~ f 
dp dz 







(29) 



Since the integral of V • (f(p(x,t), z(x,t)) u) over the computational domain vanishes thanks to the boundary 
conditions, this equation is the relation we are seeking, provided that the free energy is such that the following 
relation holds: 

df M df 
dp dz 



We are going now to reproduce this computation at the discrete level. 
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Theorem 3.1 (Stability due to the pressure work). Let C be an open convex subset of M 2 and /(•,•) be a 
convex continuously differentiable function from C to K. We suppose that (pk)kgm, {Pk)k&Mi ( z K)i<eM & n d 
{z* k )kzM are four families of real numbers such that, \/K S M., {pKi z k) £ C, (p* K , z* K ) £ C and the following 
relations hold: 



\K\ 
St 

\K\ 
St 



(pk ~ P* K )+ E Va < K P" = 

a=K\L 

(Z K - Z* K ) + ^ v t-K Za = 
a=K\L 



(30) 



where p a and z„ are given by p a — pk and z a = zk if ^a,K > 0, p a = pr, and z a = zl otherwise. Then the 
following estimate holds: 



E 

KeM 



-PK 



E 



=K\L 



> 



E i*i 



KeM 



f(p K ,Z K ) - f{P*Ki z *K) 
St 



where the family of real numbers (pk)k€M is given by: 

df df 
£ M, p K = Pk ~^-{pk,z k ) + z K -^-{pk, z k ) - f(pK,z K ) 
op oz 

Proof. Let us multiply the first relation of ((30|) by the derivative with respect to p of /(•, •), the second one by 
the derivative with respect to z of /(•, •), both being evaluated at (pk, zk), and sum: 



\K\ 
St 



( a -l\ 



(p K , 2A') V ' {PK,ZK), 



df\ \- 

OP J ( PK .z K ) a=K]L 



%) E v ^^ = ° 



(31) 



The second term of the previous relation, Tdi V ,K, can be recast as: 



df 



df) 
dz ) 



(Pk,Zk) 



^ V a ,K (Per - Pk)+ PK ^ V<T > K 
o=K\L o=K\L 

^ v ct ,k {z a - z K ) + z if E Va ' K 

a=K\L a=K\L 



(32) 



This relation is the discrete equivalent to equation (|28p : up to the multiplication by l/l-K], the first summations 
in the first term and the second term at the right hand side are the analogue of u ■ Vp and u ■ Vz respectively, 
while the second summations are the analogue of pV • u and zV • u respectively. Adding and substracting 
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f(pK, zk), we obtain a discrete equivalent of relation (|29)) : 



T, 



div.K 



= [ iL ) E y °> k (p° - p«) 



dp J 

dj_ 

dz 



{Pk.Zk) a =K\L 

^ V CT ,K (z a - Zk) + f(pK, Z K ) ^ Y a,K 
(pk,Zk) a =K\L a=K\L 



Ph \ + Z K ( ^- J - f{pK, Z K ) 



E V ^ 



In the last term, we recognize, as in the continuous setting, pk ^2 a= K\L v ^,k- The process will be completed if 
we put the first three terms of the right hand side in the divergence form. To this end, let us sum up the term 
Tdiv k over K £ M. and reorder the summation: 



^ Tdiv.K — ^ Pk 

KEM KEM 



E V "< K 
a=K\L 



E T < 



div,< 



(33) 



<jE£ h 



where, if a = K\L: 

Tdiv,a = V CT K 



d p){p K ,ZK) 



5/\ 



(Zc - zk) + /(pk, Z K ) 



(Pk,z k ) 



., J {Pc-PLj- {Za - Z L ) - f{p L , Z L ) 

In this relation, there are two possible choices for the orientation of a, i.e. K\L or L\ K; we choose this orientation 
in order to have v„,k > 0. The function (p, z) i— > f(p, z) is by assumption continuously differentiable and convex 
on the convex set C containing both (pk, zk) and {pl, zl), so the technical lemma I3TS1 hereafter applies and 
there exists (p a , z a ) in the segment [{pk , z k), (plj z l)] (itself included in C) such that: 

if (Pk, z k ) ^ (Pl, z l ) : 

) [pa - PK) + ( ) (Za ~ Z K ) + f(j)K, Z K ) 

^- [Pa ~ PL) + [^-] {Za ~ Z L ) + f{p L ,Z L ) 

otherwise: (p a , z a ) = (p K , z K ) = {pL, z L ) 

By definition, the choice (p a ,z a ) — (pcr,Zo-) is such that the term Tdiv.a vanishes, which means that the first 
three terms at the right hand side of equation (|32f are a conservative approximation of the quantity V • (/it) 
appearing in equation ()29|1 , with the following expression for the flux: 



Fa 



K 



K-. 



with: 



f & = [ ^-j {p<? - Pk) + ( %~ ) (z a - z K ) + f(p K , z K ) 

- d Pj( PK ,Z K ) \ dZ h P K,Z K ) 

.. i (p<t - Pl) + ( %- J %- z L ) + f(p L ,z L ) 



= (8J_\ 



( P L,Z L ) 
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Then the term Tdiv a can be rewritten as: 



T, 



Vo.K {Pa - Pa) 



dp) 



(PK,z K ) 



_(d_£\ 
dp) 



(pl,zl) 



With the orientation taken for er, an upwind choice yields 

Tdiv,cr = v a ,K (pk ~ Pa) 



dp) 



(PK, z K ) 



dp) 



{Pl,zl) 



+ Va,K (z K - Za) 



dz ) 



d£ 

dz 



{pk,zk) 



(Pk,z k ) 



df\ 
dz ) 



d£ 

dz 



(PL,Z L ) 



{pl-zl) 



and, by the inequality of lemma [3721 hereafter. Tdiv, a can be seen to be non-negative. Let us now turn to Tg/g t 
As the function (p, z) > f(p, z) is convex on the convex set C and both (pk, zk) and (p* K , z* K ) belong to C 
we have: 

f(PK, z K ) - f(p* Kl z* K ) 



T a/a t > \K\ 



Sf 



(35) 
□ 



Then, summing for K G M. and using relations (|3T|) . (|3"3"]) and (|33|) concludes the proof. 

In the course of the preceding proof, we used the following technical lemma. 

Lemma 3.2. Let C be an open convex subset of M. 2 , /(•,•) be a convex continuously differentiable function 
from C to K and z-{) and (p2, z?) be two distinct elements of C. Then there exists C € [0,1] such that 
(p, z) = (1 — C) (pi, Z\) + C (P2j z 2 ) satisfies the following relation: 



/(pi>*0 + (lr) 

\ d Pj(p 1 ,z 1 ) 

f(P2,Z 2 

In addition, the following inequality holds: 
T=( P i- p) 



(p-pi)+[^f) (z-zi) 



dp 



(p- pi) + 



d£\ 
dp), 



d£\ 
dp) 



(P2,Z 2 ) 



(P2,Z 2 ) 



+ {zi - z) 



dz ) 



(36) 



- Z2) 



(P2,Z 2 ) 



dZ /(puZ 1 ) \ dZ )(p 2 ,Z 2 ) 



> 



(pl, zi) 

Proof. Let us consider the function g(-) defined by: 

c~/((i-o (pu zi)+t (P2, z 2 )) 

By assumption, the function g(-) is defined over [0, 1], convex and continuously differentiable. Moreover, it may 
be checked that equation (I36| equivalently reads: 



g(0)+g'(0)t = g(l) + g'(l) (( - 1) 

or, reordering terms: 

[g'(l)-g'(0)] ( = g(0)-(g(l)-g'(l)) 
As g(-) is convex, if g'(l) = g'{0), the function g(-) is affine and g(0) — (g{l) — <?'(1)) vanishes, so the preceding 
relation is satisfied with any value of £■ Otherwise, the preceding relation allows to compute £ and, still by 
convexity of g(-), both g'(l) — g'(0) and g(0) — (g(l) — ^'(l)) is positive, and so is (. Still in this second case, 
this relation equivalently reads: 



[<?'(!) -<?'(0)] (t-l)=g(Q)+g'(0)-g(l) 
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which, as g(0) + <?'(0) — g(l) is negative, shows that ( < 1. Finally, the quantity T simply reads C [<?'(1) ~ ff'(0)L 
and is thus non-negative. □ 

Remark 3.3 (Discretization of the convective terms and conservation of the entropy). From the above compu- 
tation, it appears that the choice of p a and z a defined by equation (|34p , for the convective terms in the mass and 
the gas mass balance equations, is a convenient one to obtain an exact discrete counterpart of the continuous 
identity (|25|) . and thus, in fine, to build a scheme exactly conserving the entropy. The upwind choice yields a 
dissipation, and nothing can be said for the centered one. 

3.2. The case of a constant density liquid and an ideal gas 

Let us suppose that pe is constant and p g is linearly increasing with the pressure: 

where a is a positive real number (from a physical point of view, it is the sound velocity in a pure gaseous 
isothermal flow). For any positive p and z such that z — p + pe > 0, the relation (|6|) giving the mixture density 
as a function of the gas mass fraction and the phasic densities may be recast under the following form: 

lp = g^(p,z)= 2 91 ( 3? ) 
a z u z + pi — p 

Let us define the volumetric free energy of the mixture by: 

f(p,z) = a 2 z log(g^(p,z)) (38) 

This function is continuously differentiable over the convex subset of R 2 : 

C = {(p, z) e K 2 s.t. p > 0, z > 0, z - p + p e > 0} (39) 

We are now going to show that it verifies the other two assumptions of theorem 13. 11 namely that /(•) is convex 
and satisfies the identity: 

dp dz 

This latter relation can be proven without referring to the specific form of /(•, •), making use of the following 
property, which would be verified by the volumetric free energy function associated to any mixture composed 
of a constant density liquid phase and a barotropic gaseous phase: 

f(p,z)=zf g (g^(p,z)) with: f' g (s) = ^- 

where p(-) is the function giving the pressure as a function of the gas density (thus, in particular, f'(p g ) = p/ ' p 2 g ) 
and fg(-) stands for the specific free energy of the gaseous phase. Developping the derivatives and using the 
definition of f g (-), we get: 



T v = P — + z — ~ f = P 



dg p ' z dg p > z p 

t p = PzfgiPg)—^— + z 2 f'g(pg)—£r + z fg(Pg) ~ z f g {p g ) = z — 
op oz p 



dgP> z dg 
dp 



p,z 



+ z- 



(40) 



From the expression ([3T]) . we have: 



M1 = A_ and g§r = (41) 

dp pg z dz pe z l 
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Substituting in (|40[) leads to: 

t p = P 2 gf' g {Pa) =P 
For proving the convexity of /(•, •), we return to its explicit form: 



f(p, z) = a 2 z log 



z pe 



z + pi- p 

Differentiating twice this expression, we get: 

d 2 f 2 9 2 / 2 ( P e - P) 2 3 2 / d 2 f a p t 

— a - -7T, -—7 = a ; - - = — — — = a 



dp 2 (z + pi + p) 2 ' dz 2 z(z + pe + p) 2 ' dpdz dzdp (z + pe + p) 2 

It is thus easy to check that the determinant of the Hessian matrix A of /(•,•) is zero while its trace is positive. 
One eigenvalue of A is thus zero and the second one is positive, and /(•, •) is convex. 

4. Stability analysis 

The aim of this section is to provide some results concerning the stability of (i.e. the conservation of the 
entropy by) the scheme considered in this paper. First (section |4.1|) . in the case where both the drift velocity 
u r and the diffusion coefficient for the mass fraction of the dispersed phase D vanish (i.e. for the homogeneous 
model), we prove that the entropy (i.e. the usual entropy associated to the homogeneous model) is conserved 
by the scheme, up to a step of renormalization of the pressure which is precisely stated. Note that this step, 
which was implemented for monophasic flows in [14], could be added in the present scheme; however, we have 
chosen not to consider it further than in this theoretical section, as, in practice, its beneficial effects were not 
clear. Second (section [4. 2|) . we show that, as in the continuous case, if the drift velocity is proportional to the 
gradient of the pressure, the drift term is dissipative with respect to the same entropy; for this property to hold, 
a particular discretization of the drift term has to be implemented. 

In this section, we use the following discrete norm and semi-norm: 
VveW h , \\v\\l p = ]T \D a \ p.kl 2 

<r££int , . . 

v«ex fcl \ q \l p = £ ±^f-( qK ~ qL ) 2 

where p = (pa)ae£ in t 1S a family of positive real numbers. The function || • || 2 defines a norm over Wh, and 
| • \h, P can be seen as a weighted version of the H 1 semi-norm classical in the finite volume context [12]. The 
following relation links this latter semi-norm to the problem at hand: 

VqeL h , (BM- 1 B t q,q) = \q\l p (43) 

where B*, B and M p are the discrete gradient, (opposite of the) divergence and mass matrix defined is section 
A proof of this equality can be found in [14, section 3.4]. 



4.1. First case: u r = 0, D = 

With a zero drift velocity and a zero diffusion coefficient, the numerical scheme at hand reads, in the time 
semi-discrete setting: 
1 - solve for u n+1 

n n f.n+l _ n n-l ,.n 

" — + V • (p n u n ® u n+1 ) + Vp n - V • r(u n+1 ) = f v l+1 (44) 

ot 
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2 - solve for p n+1 , u n+1 , p n+1 and z n+1 



St 

gP' z (p n+1 , z n+1 )-p r ' 
St 

z n+l _ p n y n 



+ V(p n+1 - p n ) = 



+ V • (g p ' z {p n+1 , z n+1 ) u n+1 ) = 



(45) 



St 



V • (z n+L u n ) = 



3 - solve for y n+1 



p n+1 y n+1 = z n+1 



(46) 



Proposition 4.1 (A partial stability result). Let the density of the liquid phase be constant, the gas phase 
obeys the ideal gas law and f(p, z) be the corresponding volumetric free energy of the mixture, defined by (|38|) . 
We suppose that the viscous term is dissipative (i.e. Vw € Wh, ad(v,v) > 0). In addition, we assume that the 
density p n is positive and the gas mass fraction y n belongs to the interval (0,1]. Let u n+1 , u n+1 , p n+l , z n+1 
and p n+l be a solution to equations (I44|) - (|45[) . whith a zero forcing term. Then the following bound holds: 



~IK +1 llL« 



st 2 



+ / f(p n+ \z n+1 )dx + Sta d (u n +\u n + 1 )+"—\p n+1 \l pn 



< 



^IKIlL- 1 + Jj(p r \P n y n )dx + d ^-\p n \l 



(47) 



Proof. Multiplying each equation of the first step of the scheme (|44|) by the corresponding unknown (i.e the 
corresponding component of the velocity u n+l on the corresponding edge a) and summing over the edges and 
the components yields, by virtue of theorem 12. II 



1 

2St 



\u n+ %, p n - — +a d {u n +\u n +')- I p n W-u n+1 dx<Q 



(48) 



On the other hand, the first relation of system equation (|45|) reads, in algebraic setting: 

^M p » {u n+1 - u n+1 ) + B* {p n+1 - p n ) = 

— 1/2 

Reordering this relation and multiplying by M pn (recall that M p ™ is diagonal), we obtain: 



i M% 2 u n+1 + m;„ 1/2 b V +1 = i M% 2 u n+1 + m:„ 1/2 b V 

St St 



Squaring this relation gives: 



1 m!£V +1 + M7„ 1/2 B* p n+ \ ~ M^V +1 + M7,! /2 B* p n+1 ) = 
St 1 1 St 1 1 J 

~ m*(V +1 + m;„ 1/2 B t P n , - m^V +1 + m;„ 1/2 By 



ft 
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which reads: 



M pn u n+1 , u n+1 ) + (M-„ 1 B t p ,l+1 , B*p n+1 ) + |- BV+ 1 ) 



Multiplying by ft/2, remarking that, Vw € W/i, (M p «u, v) — ||t!||^ « and that, thanks to relation ([43]), 
V^I,,, (M~n B* q, B f q) = (BM;iB f M ) = \q\ 2 hf>ni we get: 

2 ^ 2 1 „ ' St - - . ( 49 ) 



^ ll u llh.,/3" ~~ ~2 W \h,p- 



«+l||2 " |_n|2 _ ^-.n+l ni „n 



The quantity — ({t n+1 , B*p") is nothing more than the opposite of the term / p"V • u n+1 dx appearing h 

Jn,h 

so summing (|4"5]l and (14*9")) makes these terms disappear, leading to: 



2ft 11 " ~2ft l|U ' p " 



-|b n+1 lL»-|b"lL« + (^ +1 ' Bt f n+1 )< 



Finally, (u" +1 , B*p n+1 ) is precisely the pressure work which is likely to be bounded by the time derivative of 
the volumetric free energy of the mixture. We know from theorem IA.4I that any solution to the system (|45[) 
satisfies p n+1 > 0, z n+1 > and p n+1 > 0. In view of the different forms of the equation of state gathered in 
(155]) . this implies that this solution belongs to the convex set C defined by (f3"9"]) . inside which the free energy is 
well defined, regular and convex. Hence, with this solution, theorem 13.11 indeed applies and we get: 

^ IK +1 HL>" +a d (u n+ \u n+1 ) + | \p n+ % pn + j t Jj{p^\z n+1 )d X 



< 



which concludes the proof. □ 

Theorem 4.2 (Stability of the scheme, case u r = T> = 0). Let the density of the liquid phase be constant, the 
gas phase obeys the ideal gas law and f(p, z) be the corresponding volumetric free energy of the mixture, defined 
by (|38p . We suppose that the viscous term is dissipative (i.e. Vu € Wh, a d{ v i v ) > I n addition, we assume 
that the initial density is positive and the initial gas mass fraction belongs to the interval (0,1]. 



We now add to the scheme (|44| - (|46l) the following renormalization step of the pressure, to be performed at 
the very beginning of the time step, before the velocity prediction step: 

Solve for p n+1 : -V • ( — Wp n+1 ] = -V ■ | 1 Vj>" | 

\P n J \y/p n P n ~ 1 J 

or, in algebraic setting: 

BM" 1 B* p n+1 = BM^ B* p n 

1 v p" p" 1 

Accordingly, the pressure used in the velocity prediction step must be changed to p n+1 . 
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Let (u n )o< n <N , (w™)o<n<JV , {p n )o<n<N , (^™)o<Ti<A r and (p n )o< n <N be the solution to this scheme, whith a 
zero forcing term. Then the following entropy conservation result holds for < n < N : 



\\\u n+1 \\l P ~ 



n+l 



f(p 



n+l z n-\ 



St 2 



0|2 

h,p° 



-i) dx + 5tJ2a d (u\u k )+ U -^\p^\l pn 
fc=i 

<lh°\\l P o + J^f g ( P , z °)d x + 6 -^\ P ' 

Proof. By the same proof as for the scheme without the pressure renormalization step, we get: 



(50) 



1 

2St 



^IK +1 |lL>« +a d (u n+ \u n+1 )+°i\p n+1 \l 



and the conclusion follows by summing over the time steps, remarking that z 
the renormalization step (see [14] for a detailed computation): 

\P n+1 \l P ~<\p% p ^ 



n+l 



pn+iyn+i anc ^ th a nks to 



□ 



Note that a similar pressure renormalization step has already been introduced for variable density incom- 
pressible flows [18]. 

4.2. Dissipativity of the drift term 

We address in this section the case where the drift velocity is given by the Darcy-like closure relation Q: 

1 \ Qg{p) ~ Pi r, 

u r = — (1 — a g ) a g — Vp 

A p 

In this relation, A is a positive phenomenological coefficient and a g is the void fraction, which can be expressed 
as a function of the unknowns used in the scheme as a g = z/g g (p). We recall the spatial discretization of the 
drift term in the correction step for the gas mass fraction y, namely V • (p y (1 — y) u r ), given in section[ 



G tK9(yK,yL)-G aK g(y L ,y K ) 



-K\L 



where the function g(-, •) corresponds to an approximation of <f{y) = max[y (1 — y), 0] by a monotone numerical 
flux function, K = max(C? CTj if , 0), G~ K = — mm(Go- t K, 0) and G at K is an approximation for the flux of p u r 
through the edge a = K\L. With the closure relation ([7]) for u r , a natural discretization for this quantity reads: 



Gfj.K — M Per,up 



^^(Pv-Pt) 



(pk ~Pl) 



(51) 



where pa,up is a density on a, for which, for pratical implementation reasons, we choose an upwind approximation 
with respect to the mean velocity u. The goal of this section is to show that it is possible to approximate: 



u g (1 - a g ) 



(Pg ~ Pi) 
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in such a way that this drift term is dissipative with respect to the entropy of the system. 

We begin this section by stating a consequence of the equation of state for the mixture which is central to 
the present development. 

Lemma 4.3. Let the density of the liquid phase be constant, the gas phase obeys the ideal gas law, f(p, z) be 
the corresponding volumetric free energy of the mixture, defined by (|38|) . and h(p, z) be the partial derivative of 
/(•, •) with respect to the second variable z. Then the following results hold: 

(1) h(-, ■) only depends on the pressure, i.e. there exists a function h p {-) such that, for p and z in the convex 
set C defined by (|39| . h(p, z) = h p (p(p, z)), where p(-, ■) is the function giving the pressure as a function 
of p and z: 

P = p(p, z ) = ^ 9g' Z (p, z ) = ^ — ; : 

z + pi - p 

(2) the derivative of h p (-) is given by: 

K(p) = ^# 

pePglP) 

(3) for any positive real numbers pi and pi such that pi <P2, there exists p\^ € [Pi,P2] such that: 

p ' Pi-pi 



Proof. As (p, z) E C, the pressure or, equivalently, the gas density p g can be expressed as a function of (p, z) by 
Pg = Qg' z (p, z). By the definition of /(■, •), we thus have: 



h(p,z) 



df 

dz 



fg(Pg) + z 



dz 



+ z 



df g der 
dp g dz 



Then using the expression (|4"Tj) of the derivative of qP' z (-, •) with respect to the second variable, we get: 

flip, z) — a log + p 

\a z / Pf z 



Using the fact that p = (1 — a g )pi + a g p g and thus pe — p = a g (p£ — p g ) = — (pe ~ pg), we have: 

Pg 



h(p,z) = a 2 log (J^) 



ft- Pg 
Pi Pg 



By definition of p g , i.e. p g — p/a 2 , we thus get: 



h(p, z) = a 2 



log (?) 


! p e -p/a 2 ' 






Pe 



hp{p) 



Taking the derivative of this relation yields the desired expression for h' p (-) and, as h p (-) is continuously diffcr- 
entiable in [pi,P2], the existence of p\^ follows by Lagrange's theorem. □ 



We are now in position to state and prove the following stability result. 
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Proposition 4.4. Let the density of the liquid phase be constant, the gas phase obeys the ideal gas law and 
f(p,z) be the corresponding volumetric free energy of the mixture, defined by (|38[) . Let C be the convex set 
defined by ([55]) and (pk)keM> {"}Jk)k£M an d { z k)k^m be such that, \/K £ M, (pk,PkUk) S C, (pk,zk) G C, 
and the following relation is satisfied: 



-jjj- (pk Vk - z K ) + G+ K g(y K ,yL) - G~ K g(y L ,yK) =0 

a=K\L 



(52) 



where g{-, •) corresponds to an approximation oftp(y) = max[ y (1—y), 0] by a monotone numerical flux function, 
K = max(Gcr t x, 0), G~ K = — min(G CTi K, 0) and G a ^K is given by the relation (jSTj) . Then, if g(yK,VL) > 
for all a £ £i n t, o~ = K\L, there exists a discretization for the term: 

a g (1 - a g ) . 
- JL —, (Pg - Pi) 



in (|5ip such that the following stability estimate holds: 

I \ K \ if(pK,PKyK)-f(pK,Z K )}<0 



KeM 



which means that the drift term is dissipative with respect to the entropy of the system. 

Proof. We multiply equation (|52"|) by the partial derivative of /(•, •) with respect to the second variable, taken 
at the point (pk, Pk Vk), and sum up over the control volumes of the mesh: 



Y Kpk,PkVk) 



KeM 

where T\ and Ti reads: 



\K\ 

St 



(PkVk ~ %k)+ Yl G a,K9iyK-,VL) - G a K g(y Ll y K ) 



=K\L 



Ti + T 2 = 



\K 



Ti = Y Kpk,PkVk) [pkVk-zk] 



KeM 



T 2= Y KPK,PKVk) 

KeM 



G+ K g(y K ,yL)~G aJ( g(y L ,yK) 

a=K\L 



As the fonction /(•,•) is convex, we have: 

Ti > t- Y \ K \ \I(Pk,PkVk) - f(PK,z K )} 



(53) 



KeM 



Let us turn to T%. Reordering the sum, we get 

Tz= } J \<y\ Pa,u P gu P (yK,yL^ 



tre£i, 



Ug (1 - OLg) 



(Pg - Pi) 



(PK-Pl) [Hpk,Pk Vk) - Kpl-.Pl Vl)\ 



where g xl . v {yK,VL-,u r ) = g{yK,VL) if u r > and g np {VK,VL,u r ) = g{yL,VK) if «r < 0; in any case, we have, 
by assumption, g U p(yK,yL,Ur) > 0. We now choose, for the approximation of the quantity defined on a in the 



2G 
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preceding relation, an expression of the form: 



[Pg - Pi) 



A 



A 



where (a g ) a stands for an approximation of the void fraction on a which only needs here to be supposed 
non- negative. Applying lemma [4731 T 2 reads: 



T 2 = ^2 W Pa ' u P 9n P {VK,VL,u r ) <ya ^ a — p e Qg{p a ) h' p {p a ) (p K ~ Pl) [h p {p K ) - h p (p L )} 

If Pk — Pl, the term associated to K\L in this sum vanishes. Otherwise, from the third assertion of lemma 
14.31 there exists p a € [mm(p Kl PL), max(p x , p L )] such that the product h' p (p a ) (pk — Pl) [h p {pK) ~ h p (p L )] is 
positive. Since we choose (a g ) a such that {a g )a > 0, all the other quantities are positive, and this concludes 
the proof. □ 

The following proposition extends the stability result of the preceding section to the case u r 7^ 0. 

Proposition 4.5 (Stability of the scheme, case u r 7^ 0). Let the density of the liquid phase be constant, the 
gas phase obeys the ideal gas law and f(p, z) be the corresponding volumetric free energy of the mixture, defined 
by ()38p . We suppose that the viscous term is dissipative (i.e. Vi> € Wh, ad(v,v) > 0). In addition, we assume 
that the density p n is positive and the gas mass fraction y n belongs to the interval (0, 1]. Let u n+l , u n+l , p n+l , 
z n+ , p n+1 and y n+l be a solution to the equations of one time step of the scheme, whith a zero forcing term. 
We suppose that the drift velocity is given by a the Darcy-like relation ([7]) and that the discretization of the 
correction step for the gas mass fraction y n+1 is such that the stability result of proposition \4-4\ applies. Then 
the following bound holds: 



1 r st 2 

1 



^ h n+1 \\i p n + I f(p n+ \p n+1 y n+1 ) dx + Sta d (u n+1 ,u n+1 ) + —\p n+1 \l pn 



< 5 h n \\i P « 



r /)t 2 
1 + Jj{p n ,p n y n )dx + —\p n \l. pn 



Proof. Proposition 14.41 yields: 

i E 1^1 [f(PK + \p n K +1 y n K +1 )- ftP n K^ n K )] <0 

KeM 

The conclusion thus follows by summing this relation with the estimate of proposition 14.11 □ 

Finally, note that, as in the preceding section, this partial stability result yields the same entropy decrease 
estimate for the whole scheme as in the preceding section if a renormalization step for the pressure is added to 
the scheme. 

Remark 4.6 (On the choice of the monotone numerical flux function). As stated in section[2l we have adopted 
for the numerical tests presented hereafter the following flux-splitting formula: 

g(a 1 ,a 2 ) = ffi(ai) + .9 2 (a 2 ) 

where 31 (ai) = ai if ai £ [0, 1] and zero otherwise and and 32(02) = — (a 2 ) 2 if «2 £ [0, 1] and zero otherwise. 
This numerical monotone flux does not satisfy the hypothesis of proposition ^. 51 as it is not always non-negative. 
However, several other choices are possible for the numerical flux function g(-, ■) (e.g. [12]), and some of them 
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solve this problem. Thanks to the fact that ip(s) = s (1 — s) is positive Vs € [0, 1], it is the case, for example, 
for the flux obtained with a one-dimensional Godunov scheme for each interface: 



0(01,02) 



max{</?(s), 0,2 < s < a±} if 0,2 < a,i 
min{<yj(s), ai < s < 02} if ai < 02 

5. Numerical results 



This section is devoted to numerical tests of the proposed scheme. We first adress a problem built in such 
a way that it admits an analytical solution, to assess the convergence properties of the scheme. Then several 
additional tests are performed, to check the stability of the algorithm and the quality of the results. 

5.1. Assessing the convergence against an analytic solution 

We address here a problem built by the so-called technique of manufactured solutions: the computational 
domain and the solution are chosen a priori and the initial conditions, the boundary conditions and the forcing 
terms are adjusted consequently. Let thus the computational domain be = (0,1) x (—1/2,1/2), and the 
density and the momentum take the following expressions: 



p = 1 + — sin(7rf) [cos(7ra;i) — sin(7ra2)] pu = — — cos(nt) 



sm(7ra;i) 

C0S(7TX2) 



The pressure and the partial gas density are linked to the density by the equation of state (fTT|) , where the liquid 
density pe is set at pi = 5 and the quantity a 2 in the equation of state of the gas (O is given by a 2 = 1 (so 
p g = p). We choose the following expression for the unknowns y and z: 

2.5- 0.5 p 2.5- 0.5 p 

y = — 7~5 z = py= — 

4.5 p 4.5 

The relative velocity is constant and given by u r = (0, 1)* and the diffusive coefficient D is set to D — 0.1. The 
analytical expression for the pressure is obtained from the equation of state (i.e. relation ([37])). These functions 
satisfy the mass balance equation; for the gas mass fraction and momentum balance, we add the corresponding 
right-hand side. In this latter equation, we suppose that the divergence of the stress tensor is given by: 

V't(u)=/iAu+^VV'W, p = 10~ 2 

and we use for the viscous term the corresponding form for the bilinear form ad(-, •) (see section [23]) • 

Errors for the velocity, pressure and gas mass fraction obtained at t = 0.5, as a function of the time step and 
for various meshings, are drawn on figure [2] figure [3] and figure [U respectively. These errors are evaluated in the 
L 2 norm for the velocity and in the discrete L 2 norms for the pressure and the gas mass fraction. Computations 
are made with 20 x 20, 40 x 40 and 80 x 80 uniform meshes (so with square cells and the Rannacher-Turek 
element). For large time steps, these curves show a decrease corresponding to approximately a first order 
convergence in time, until a plateau is reached, due to the fact that errors are bounded by below by the residual 
spatial discretization error. The value of the errors on this plateau then show a spatial convergence order close 
to one, which is consistent with the choice of an upwind discretization for the advection terms in the mass and 
gas mass fraction balance equations. 

5.2. Two-dimensional sloshing in cavity 

Two layers of non-miscible fluids (air and water) are superimposed with the lighter one on top of the heavier 
one. The gravity (with g — 9.81 m.s~ 2 ) is acting in the vertical downward direction. The length of the 
rectangular cavity is L = lm, the height of each layer is respectively hi = lm and h g = 1.25m, so the total 
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l.e-4 



x x Mesh 20x20 
x x Mesh 40x40 
v— v Mesh 80x80 



l.e-3 l.e-2 
Time step 



l.e-1 



Figure 2. Error for the velocity at t = 0.5, as a function of the time step (L 2 norm) 




Mesh 20x20 
Mesh 40x40 
Mesh 80x80 



l.e-3 l.e-2 
Time step 



l.e-1 



Figure 3. Error for the pressure at t = 0.5, as a function of the time step (discrete L 2 norm). 



height of the box is 2.25 m. The water and air densities are respectively pi = 1000 kg.m~ 3 and p g = p/a 2 where 
a 2 is such that p g = 1.2 kg.m^ 3 at p = 10 5 Pa. The diffusion coefficient D and the drift velocity are set to zero. 
A perfect slip condition is imposed on the whole boundary. At initial time, both fluids are at rest, then the 
cavity is submitted to an horizontal acceleration given by clq = 0.1 m.s~ 2 . 



In the case where both fluids are supposed incompressible and the convection and diffusion terms may be 
neglected, an analytical solution for the flow in a rectangular cavity is provided in [5]. In particular, the shape 
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FIGURE 4. Error for the gas mass fraction at t = 0.5, as a function of the time step (discrete 
L 2 norm). 



of the interface is given by the following relation: 



9 



L 

x 

2 



4^ Lk 



cos(w 2 „+i t) cos(fc 2 „+i t) 



2n+l 



where the wave number k n is defined by: 

2-kti 



and w„ is given by: 

^2 = 9 k n {Pi ~ Pg) 

Pg coth(fc„ h g ) + pi coth(fc n he) 
In practice, to compute this analytical solution, we perform the summation up to n — 200. 

As, to remain in the domain of validity of the solution, the amplitude of the fluid oscillations must be very 
small, a very fine mesh is necessary near the free surface, to capture its motion. The mesh is thus made of about 
41 000 rectangular cells (with the Rannacher-Turek element) and, in the vertical direction, the space step is 
adapted in such a way that it is smaller near the interface between the two phases and equal to 8x2 = 0.0005 m, 
and increases when moving away the free surface, up to 5x2 — 0.05 m at the top and bottom sections. In the 
horizontal direction, the mesh is uniform with step size 5x\ = 1/70 m. Calculations with different viscosities 
have been performed, these latter being supposed to vary with the mixture density: fi = p/100, fi = p/1000, 
fi = p/10000. 

The numerical results are reported on figure[5](^ = p/100), figure[S](^ = p/1000), and figure[7](p = p/10000) 
respectively. Comparing the obtained shape for the interface with the analytical solution, we observe that the 
numerical solution is closer to the analytical one with fi = p/1000 than with /i = p/100, certainly because 
the fluid is too viscous in this latter case. More surprisingly, when reducing the viscosity to p = p/10000, the 
numerical solution also becomes less accurate. Our explanation is that, to obtain a good solution, it is necessary 
to respect a balance between approaching the physical problem (which, in this case, would suggest p = 0) and 
keeping sufficient coercivity to ensure a reasonable convergence of the numerical approximation (which, on the 
contrary, requires a high value for the viscosity). With a more refined mesh, viscosity thus probably could be 
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decreased, and the solution be closer to the analytical one. However, with this mesh already, results seem to be 
rather more accurate as those available in the litterature [5] . 




0.5 1. 1.5 

Time 



Figure 5. Sloshing in cavity: analytical solution and numerical solution with fi = p/100. 
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Figure 6. Sloshing in cavity: analytical solution and numerical solution with /i = p/1000. 
5.3. Bubble column 

We address in this section a classical benchmark for diphasic flow solvers, namely the flow in a pseudo two 
dimensional bubble column investigated experimentally by Becker et al. [3]. The apparatus has a rectangular 
cross section with the following dimensions: its width is L = 50 cm, its depth is 8 cm and it is H — 200 cm high 
(see figure [5]). ft is filled with water up to the height h = 150 cm. A gas sparger, positioned 15 cm from the 
left wall, is used to introduce an air flow of q = 8l.min into the system. The circular sparger has a diameter 
of 40 mm and a pore size of 40 /xm. Several liquid circulation cells can be observed in the column, the location 
and size of which continously change. The bubble swarm is influenced by these vortices and therefore rises in a 
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Figure 7. Sloshing in cavity: analytical solution and numerical solution with p = pj 10000. 



meander-like way. The direction of its lower part is stable and directed towards the nearest sidewall; its upper 
part changes its shape and location in a quasiperiodic way, according to transient liquid circulations [27]. 

To simulate this experiment, we choose the following data. The boundary conditions are defined at the inlet 
as follows: 

_ q 

*^ <^(/,imp 

where S is the gas inlet area and a s ,imp = 1 is the void fraction imposed at the inlet. Along the walls and at the 
outlet of the column, homogeneous Dirichlet conditions are used for the velocity. Initial conditions are set to u = 
Oro.s -1 and p — po where po = 10 5 Pa is the ambiant pressure. The density of the liquid is pi — 1000 kg.m~ 3 ; 
the gas obeys an ideal gas equation of state p g — p/a 2 , where a 2 is such that p g = 1.2 kg.m~ 3 at p = 10 5 Pa. 
The diffusion coefficient D is set to zero, the drift velocity is constant and given by u r = (0, 0.2)* m.s^ 1 . 

For this test case, we use a regular meshing composed of rectangular cells (with the Rannacher-Turek ele- 
ment) with 76 meshes in the horizontal direction, from which 4 for the gas inlet, and 300 in the vertical one. 
Calculations with time steps up to St = 10 _1 s have been performed, observing that smaller time steps yield a 
thinner free surface. 

The viscosity is a parameter difficult to adjust, since, in this simulation which is based on the system of 
equations governing a laminar flow, it must represent in some way the turbulent diffusion, i.e. the effects 
of fluctuations of the flow at microscopic scales, which may originate from the usual turbulence phenomena 
(sometimes termed "monophasic turbulence") and from the perturbation of the velocity field due to the motion 
of the bubbles (sometimes termed "diphasic turbulence"). Calculations with a viscosity ranging from p = 
10 -3 Pa.s to p = 10 2 Pa.s have been performed. With smaller viscosities, we observe more oscillations of the 
free surface, the bubble swarm reaches the free surface faster and is farther from the sidewall. 

Finally, the numerical results obtained with St = 10~ 2 s and a viscosity of p, = 1 Pa.s are reported on figured! 
With this value of the viscosity and these mesh and time steps, numerical convergence seems to be reached, at 
least visually. One can observe the stability and the thinness of the free surface. Results qualitatively reproduce 
the expected behaviour, which is the best we can hope with the rather crude modelling of turbulence which we 
adopted. 
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H 




Figure 8. Bubble column: geometry of the problem and void fraction at times 2s, 4s and 40s. 

6. Conclusion 

In this paper, we adress the drift-flux model, which, for isothermal flows, consists in a system of three balance 
equations, namely the overall mass, the gas mass and the momentum balance complemented by an equation of 
state and a phenomenologic relation for the drift velocity. 

For this problem, we develop a pressure correction scheme combining finite element and finite volume dis- 
cretizations, which enjoys the following properties. First, the existence of a solution to each step of the algorithm 
is proven. Then, essential stability features of the continuous problem still hold at the discrete level: the un- 
knowns are kept within their physical bounds (in particular, the gas mass fraction remains in the [0, 1] interval); 
in the homogeneous case {i.e. when the drift velocity vanishes), the discrete entropy of the system decreases; 
in addition, when using for the drift velocity the Darcy-like relation suggested in [20], the drift term becomes 
dissipative. Since, when the density is constant, this fractional step algorithm degenerates to an usual incre- 
mental projection method based on an inf-sup stable approximation, stability can be expected in the zero Mach 
number limit. Finally, the present algorithm preserves a constant pressure and a constant velocity through 
moving interfaces between phases (i.e. contact discontinuities of the underlying hyperbolic system). To achieve 
this latter goal, the key ingredient is to couple the mass balance and the transport terms of the gas mass balance 
in an original pressure correction step. 

We chose in this paper to only consider the case of a constant density liquid phase and of a gaseous phase 
obeying the ideal gas law. Dealing with a more general barotropic gas phase is certainly the simplest generaliza- 
tion, but the present theory also seems to extend to the case of a compressible fluid with minor modifications: 
for the stability study, essentially, the expression for the volumetric free energy of the mixture should be replaced 
by the usual expression applying when both phases are compressible, see for instance [20] ; the existence theory 
would probably be simpler, since an upper bound for the density would provide in this case an estimate for 
the pressure. Returning to the case of an incompressible fluid, extending the present theory to deal with pure 
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liquid zones appears on the contrary to be a difficult task, since the role played by the pressure in such a system 
seems to deserve some clarifications. 

Numerical tests show a near-first-order convergence in space and time, consistent with the implemented 
discretization: first order backward Euler method in time and standard upwinding of the convection terms in 
the mass and gas mass fraction balance equations. With respect to this latter point, using more accurate space 
discretization (typically, MUSCL-like techniques) should certainly be desirable. 

To assess the robustness of this algorithm, various numerical tests have been performed. They show in par- 
ticular that free surface flows are computed without any instability, keeping a rather sharp interface throughout 
the computation. In addition, pure monophasic liquid zones are supported, although, as already mentioned, 
this case remains beyond the scope of the theory developped here. This scheme is now implemented in the ISIS 
code developped at IRSN and daily used for industrial applications. 

Appendix A. Existence of a solution to a class of discrete diphasic problems 

We address in this section the following abstract discrete problem: 

a(u, <p$) - [ p V ■ <p®te = [ f ■ <p®dx, Va e £ iut , l<i<d 

Ja,h Jn 

lQ[ e *>(p Kl z K )-Q p ' z (p* K ,z* K )]+ ]T v+ iK q**(pk,zk)-v- k Q*'(pl,z 1 .)=0, VK e M (m) 

a=K\L 

\K\ 

LJ.( ZK - Z * K )+ v+ k z k -v- k z l = Q, VK e M 

a=K\L 

This problem is supposed to be obtained from (part of) a continuous problem by a space discretization combining 
Rannacher-Turek or Crouzeix-Raviart finite elements and finite volumes; notations related to discrete quantities 
are given in section [2.21 and are not recalled here. The bilinear form &(•,•) is just assumed to be such that 
||u|| a = [a(u, u)] 1 / 2 defines a norm over the discrete space W h . The quantities (v+^)™ +1 and (v~ K ) n+1 stands 
respectively for max(v™ H ^ L , 0) and — min(v™ H ^ L , 0) with v™^ = |er| • ukl- Note that, in the last two 
equations, the flux summation excludes the external edges, which implicitely expresses the fact that the velocity 
is supposed to vanish on the boundary. 

This system must be completed by three equations of state. The first two ones giving the liquid density pi 
and the gas density p g as a function of the pressure: we suppose here that the density of the liquid is constant 
and that the gas obeys the equation of state of ideal gases, which, for the sake of conciseness, we suppose here 
to be simply p g — p. The last equation relates the mixture density p with the gas mass fraction y or the gas 
partial density z = py and the phases density, and may take the three following forms: 

P=p(p,z)= * Pl P = Q p '*(p,z)=pi + z(l-^) p = g P>y(p,y) = - I (55) 

z + p e - p p t-y y 

pi p 

These three relations are equivalent as soon as the following assumptions for the unknowns of this system are 
satisfied: 

p > 0, p > 0, z > and < y < 1 (56) 

These assumptions are natural, excepted the hypothesis that y or z does not vanish, which excludes the existence 
of purely liquid zones. This latter assumption is assumed to hold for the initial quantities, i.e. we suppose that: 

VKeM, y* K = ^e(0,l] (57) 
Pk 
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where p* K = g p ' z (p* K , z* K ). 

Our aim in this section is to prove that there exists a solution to system (1541) complemented with one of the 
relations of (|55p . under the assumption (|57p . and that any such solution satisfies the inequalities (|56l) . 

We begin this section by two preliminary lemmas. 

Lemma A.l. Let (x k )keM an d (xk)kgM be two families of real numbers satisfying the following set of 
equations: 

\K\ 

VK e M, -^-(x K -x* K )+ 22 v t,K x k - v~ K x L = 

a=K\L 

We suppose that, WK £ M, x* K > 0. Let ||V^ • it||oo be defined by: 



|Vfc-«|| 



max 

KeM 



K 



a=K\L 



Then, \/K £ M, Xk satisfies: 



mm x K 



<XK< . > , \K\ X* K 



1 + St ||Vh ■ u\\oo ~ min \K\ A-^ 

KeM KeM 



Proof. The first inequality follows from an application of the discrete maximum principle lemma which can 
be found in [14] (lemma 2.5, section 2.3). The second one then follows from the fact that, by conservativity, 
J^KeM Xk ~ J2k£M x K> remarking that, by the preceding relation, the values xk, for K £ A4, are all 
positive. □ 

The proof of the following result can be found in [23] . 

Lemma A. 2. Let {p* K )KeM, (x* K )i<eM, {pk)k&M an d { x k)k^M be four families of real numbers satisfying 
the following set of equations: 

\K\ 

VKeM, — {p K X K - p* K X* K ) + V i. K PK X K - V~ K PL X L = 

(7—K\L 

We suppose that, \/K £ M, p* K > 0, px > and: 

o=K\L 

Then the following discrete maximum principle holds: 

WK £ M., min x* T < xk < max x* T 
LeM LeM L 



We now state the abstract theorem which will be used hereafter; this result follows from standard arguments 
of the topological degree theory (see [8] for an overview of the theory and e.g. [11,14] for other uses in the same 
objective as here, namely the proof of existence of a solution to a numerical scheme). 
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Theorem A. 3 (A result from the topological degree theory). Let N and M be two positive integers and V be 
defined as follows: 



v = {(i 1 j,z)ei iV x 



such that y > and z > 0} 



where, for any real number c and vector y, the notation y > c means that each component of y is greater than 



c. Let b £ 



pAI 



to M JV x R x R satisfying: 



and /(•) and F(-, ■) be two continuous functions respectively from V and V x [0, 1] 



(i) = /(•); 

(ii) V# € [0, 1], if an element v of O (the closure of O) is such that F(v, 9) = b, then v G O, where O is 
defined as follows: 



O = {{x,y,z) £ R iv x 



s.t. \\x\\ < M and e <y < M and e < z < M} 



with M and e two positive constants and || • | a norm defined over M. N ; 
(iii) the topological degree of F(-,0) with respect to b and O is equal to rf 7^ 0. 

Then the topological degree of F{-, 1) with respect to b and O is also equal to do ^ 0; consequently, there exists 
at least a solution v € O such that f(v) = b. 

We are now in position to prove the existence of a solution to the considered discrete system. 

Theorem A. 4 (Existence of a solution). Under the assumption (|57[) , the nonlinear system (|54[) complemented 
with the relation (155} admits at least one solution, and any possible solution is such that: 



VK&M, 



p K > 0, z K > 0, < y K 



zk 

PK 



< 1, 



p K > 



Proof. This proof makes use of theorem IA.3I twice, by linking the initial problem to a linear one through two 
successive homotopies. Let N — d card(£j nt ) and M = card(A / I); we identify the finite element space of discrete 
velocity with M. N and the finite volume space of pressure and partial density with R M . Let V be defined by 
V = {{u,p, z)eR N x R M x R M such that p > and z > 0}. 



Step 1: first homotopy 

We consider the function F : V x [0,1] 



given by: 



F(u,p,z,6) 



sk 



a(u, ¥>«)-/ p V • ¥>W Ax- f- vjW dx, 
Jn,h Jn 

\K\ 

-^r [q!' z (pk, zk) - Qe' z {p K ^ z *k)\ 



er G £ int , 1 < i < d 



E 



(jPK,ZK) -V K Qf (PL,Z L ), 



= -^f[ZK ~ Qe ,Z (PK^ z K)V*K\+ E Y t,K Z K -V a K Z L , 



-K\L 



St 



a=K\L 



KeM 



(58) 



where the function Q@'' z (■, ■) is implicitely defined by the following relation: 



i - y , y 

Qifiip) P 



with ge,e(p) 



1-9 



and z — py 
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Note that this definition makes sense (i.e. using z — py, the function g p ' z (-, •) can be explicitely computed from 
the expression of g p ' v (-, •)) as soon as p > 0, and thus for any (u,p, z) £ V. 

Problem F(u,p, z, 1) = is exactly the same as system (|54p. 

Let e and M be two positive real numbers, and O be defined by: 

O = {(u,p, z) e R N x R M x K M s.t. ||u|| a < M, e < p < M and e < z < M} 

We now suppose that (u,p,z) S O (and thus, in particular, p > e) and that F(u,p, z,9) — and provide 
estimates for (u,p, z). 

We begin by the following elementary bound, which is useful throughout the proof. From the definition 
°f 8e,e(p), we observe that mm(pe,p) < gi.e(p) < m.8tx(pt,p). By the same way, provided that y £ [0,1], 
mm(Qg t e(p),p) < gg' z (p,z) < max(ge^(p),p). Hence, mm(pe,p) < gg' z (p,z) < max(pe,p) and, thanks to 
assumption ([57]) : 

V9 e [0, 1], VK € M, ^"'(px, 4r) < P* with p* = max 



(max^), Pi 



Step 1.1: || • || a estimate for the velocity. 

Let us first recast the equation of state of the mixture under a more convenient form. Substituting its definition 
for g e ,e(p) in ge' v (p,y), we get: 

9 = 9(l-y) | y + (l-8)(l-y) = 1 - y> | g (59) 
Qi P Qi P 

with y'(y, 9) = y+(l—6)(l—y). Then, taking y — zj gg ,z {p, z) as unknown in the third equation of F(u,p, z, 9) 
0. we get, for any K G M: 

-jjT [Qe° ,Z (PK,z K ) y K - gg°' z (p* Kl z* K ) y* K ] + ^ v+ K g^ z (p Kl z K ) y K - v~ K gg ,z {pL,z L ) y L = 

<j=K\L 

As, by the second equation of F(u,p, z, 9) — 0, this relation vanishes for the constant function yx — 1, V-ftT G M, 
we also obtain: 

-JjT [q$ ' z (pk, z k ) y' K - gf z (p* K ,z* K ) (y')* K \ + ^ v+ K g^ z (p K , z K ) y' K ~ v~ K gf z {pL,z L ) y' L =0 (60) 

a=K\L 

where, by assumption {ST}, VK G M, (y')* K = Vk + i 1 ~ ^X 1 - Vk) e i 1 ~ 6 + 1]. with y* = mm KeM y* K . 
We thus obtain a new problem, which keeps the structure of system (l54l) . with the same equation of state (i.e. 
relation ([59]) ) and just a modified initial value for z (i.e. z* K changed to gg ,z (p^, Z * K ) (y')* K ). The unknown p is 
still an unknown of this new problem, and we thus have by assumption p > 0. In addition, by lemma lA.2[ any 
solution of this new problem is such that the gas mass fraction verifies 1 — 8 + #y* < y < 1, and thus the density 
and the gas partial density are positive. The unknowns thus belong to the domain where the free energy is 
correctly defined, and theorem 13 . 1 1 applies . Multiplying the first equation of F(u,p, z,9) = by u CT ^, summing 
over a G £j nt , 1 < i < d and using Young's inequality thus yields: 

l II II 2 l ^ ^ I ■*■ "\ n^./ \ I 1 / \ -l X ^ l-r^l «. K / Sli 4t \ / f\ * 1 / * \ l|lflll9 



Ml+i E 1^1 erkPK,z K )y' K logfc*) < \ E 1*1 Q*'Vk,*kWYk *>&(Pk) + 2 
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where || • ||_ a stands for the dual norm of || • || a with respect to the L 2 inner product. The summation at the 
right hand side of this relation is bounded by (1/St) \ p* log(p*) where p* = m&stKeM V*k- By conservativity 
of equation ([501) , J2kgm \^\ 6s(p*k^ z *k) (v'Tk — 1^1 P* ■ Since, by assumption, p > e, we thus get: 

||u||2 < | |0| p* |log(e)| + | |0| p*log(p*) + ||/|| 2 _ a 
For e > small enough, we thus have: 

HI. <ci |log(e)| 1/2 (61) 
where, in this relation and throughout the proof, we denote by Cj a real number only depending on the data of 
the problem, i.e. f2, p* , p* , /, a(-, •), St and the mesh, and the expression "e small enough" stands for e < c[ 
where c[ is a positive real number itself only depending on the data. 

Step 1.2: L°° estimates for z. 

By equivalence of the norms over finite dimensional spaces, inequality ()61[) also yields a bound for u in the L°° 
norm and, finally, for ||V/, • w||oo : 

||V/,.u|U <c 2 llog(e)! 1 / 2 

By lemma [ATTT we thus get from the third relation of the system F(u,p, z, 9) — 0, still for e small enough: 

^>c 3 |log(e)|- 1 / 2 (62) 
On the other hand, we get from the same relation by conservativity: 

z < c 4 (63) 

Step 1.3: L°° estimates for p. 

From the first relation of (j55]) . using the bounds for z, we get: 

p>c 5 llog^r 1 / 2 (64) 

To obtain an upper bound for p, we first remark that, as the considered spatial discretization satifies a discrete 
inf-sup condition, a bound for u provides a bound for p — m{p) where m{p) stands for the mean value of p. By 
equivalence of norms on finite dimensional spaces, we can choose to express this bound in the seminorm defined 
by \/q € Lh, \q\\ y \.h — J2cr££ iIIt (<j=k\l) \ok ~ Ql\- With this semi-norm, the mean value of p disappears, and we 
get for e small enough: 

J2 Ipk-PlI <c 6 |log(e)| 1/2 (65) 

<xe£ int (<j=K\L) 

An upper bound for p in one cell of the mesh, say Kq, would then provide an upper bound for p, since, for any 
K 6 M, it is possible to build a path from Kq to K crossing each internal edge at most once. To obtain such 
an estimate, we follow the following idea. If the pressure is somewhere lower than pj>, we are done; otherwise, 
with the chosen equation of state g@ :Z (-, •), when 9 varies, the liquid is everywhere denser than for 9 = 1 and we 
are going to show that, even if its total mass also increases, the volume that it occupies is lower than for 9 = 1. 
Hence, the remaining volume for the gas is bounded away from zero, and, by conservation of the gas mass, the 
pressure cannot blow up everywhere. First, we need to introduce the phase volumetric fractions. The equation 
of state ((55)) can be written as: 

P , P-z 



1 



V Pi 
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and, as p — z = p (1 — y) and y < 1, both fractions at the left hand side of this relation are non-negative. We 
may thus define a g € [0, 1] and ai 6 [0, 1], referred to as the gas and liquid volume fraction respectively, by: 

P P-z 

Oig = ~ Ott = 

V Pi 

Note that a g + at = 1. Combining the second and the third relation of the system F(u,p, z, 9) = 0, summing 
over the control volumes of the mesh and remarking that the fluxes cancel by conservativity, we get: 

E \K\ ( ae ) K = ]T \K\ g« (66) 

KeM KeM Pe 

Let us denote by (a^)K.i the liquid void fraction with the equation of state corresponding to 9 = 1: 

1-24 



v*k , 1 ~y* K 



VPk pi 

Exploiting the form (|5^|) of the equation of state for 9 ^ 0, we obtain from relation ([5b]) : 

Vk , l - y* K 

E \K\ (a t ) K = E \ K M)^ J)* K i-\y'y K with (v')k=Vk+0Q--Vk) 



KeM KeM 



Vk Pt 



If we suppose that pk > Pi-, the fraction in the above equation is bounded by 1: indeed, both the numerator 
and the denominator are harmonic averages of p* K and pi, the weight associated to p* K being larger in the 
denominator, since {y')* K is closer to 1 than y* K . We thus get: 

E 1*1 Kk>c7 = M- E l*K^W 

KeM KeM 

where C7 is positive by assumption, since VK E Ai, y* K > 0. Thus there exists Kq £ Ai such that (a g )K > 
c$ = C7/|Q| . On the other hand, we have, still by conservativity: 

E 1*1 K)* vk = E \ K \ z * = E 1*1 z *= £ \ K \ qV\v k ^k)vk < m p* 

KeM KeM KeM KeM 

We thus get, since all the {ct g )K and pk are non-negative: 

(a g ) K „ PK Q < P* 

and thus, as {ct g )K is bounded by below, the pressure is bounded by a quantity only depending on the data. 
As a consequence, for e small enough: 

p<c 9 llog(e)! 1 / 2 (67) 
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Step 2: second homotopy 

We consider the function F : Fx [0,1] 



given by: 



F(u,p,z,9) = 



V a>i = a(u,(p ( *>) 

\K\. * s 



6 f p V • <p® dx - [ f ■ cp¥ dx, 
Jn.h Jn 



SK 



\K\ 
St 



a=K\L 



P K a=K\L 



V G £intj 1 < i < d 

KeM 



(68) 



The system F(u,p, Z, 1) = is the same as the system obtained at the end of the preceding homotopy for 6 = 0, 
and the system F(u,p, z, 0) = is linear and clearly regular (by stability of the bilinear form a(-, •)). 

In addition, the third equation is now decoupled from the first two ones, and these latter have the structure 
of a monophasic compressible problem as studied in [14]. From this theory, an estimate similar to the first one 
in the preceding step is available and reads: 

\ Ml + i E 1*1 P* lo s(PA') < Yi E I*' Pk log&&) + \ \\f\W 



KeM keM 
Since the function s i— > slog(s) is bounded by below on (0, +oo), this latter relation yields: 

|M| a < Cio 

By lemma POI we thus directly get: 

Cll < p < C12, Ci 3 < Z < CM 



(69) 



(70) 



Conclusion 

We choose e small enough for the relations (j6T|) . (|6"2"]) . ((65|l and (IbTjl to hold, e < min(cn, C13) and, in addition: 

e < max(c 3 , c 5 ) I log(e)|~ 1/2 

which is possible because the function s 1— ► slogs tends to zero when s tends to zero. Let now M be such that: 



M > max 



max(ci,c 9 ) I log(e)| 1/2 , c 4 , ci , c i2 , c u 



Then, from inequalities t[6Tj). ([62]) . ([65]). (|64| . l[67]l . l[69|l and ([70]). we get that throughout both homotopies, the 
unknown (ui,p, z) remains in O. As the last linear system is regular and admits a solution in 0, the topological 
degree of F(-, •, •, 6) with respect to O and zero remains different of zero all along both homotopies, which proves 
the existence of a solution in O. 

We now turn to the proof of the a priori estimates p > 0, z > 0, < y* < 1 and p > 0. The fact that, if 
p* > and z* > 0, then p > and z > is a direct consequence of lemma POl applied to the second and third 
relation of problem ([54]). In addition, as both p* > and p > 0, lemma IA~2l applies and thus, as < y* < 1, we 
have < y < 1. If p > p^, the fact that p > is evident. In the other case, by the equation of state written as 
a function of p and y (third form of (|55|). we get first that: 



P< P < Pi 



(71) 
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and, second, that, since p > 0, the pressure does not vanish. The second form of this same relation (|55| thus 
can be written: 

z z 

P = - P + (1 - -) Pi = Otg P + (1 - OLg) Pi 

As p ^ pi, the void fraction a g thus reads: 

P~ Pi 

OLg = 

Pi -P 

which, by inequalities (|7ip . yields a g > and, finally, since z > 0, p > 0. □ 

This existence result applies directly to the pressure correction step used in the algorithm presented in 
this paper, with a particular expression for the bilinear form a(-,-), which reads, dropping for short the time 
exponents: 

a[U,V)= 2^ —<rj- Pa Ua ■ V a 

Note that the analysis is performed here with a very simple equation of state for the gas (p = p), but would be 
readily extended to general barotropic laws p — p(p), under the mild assumptions that the corresponding free 
energy exists and is convex and the function p(-) is increasing and one to one from (0, +oo) to (0, +oo). 

Let us now turn to the discretization of a stationary diphasic problem. As happens in the monophasic 
case, [16], it is likely that, in the case where the velocity is prescribed on the whole boundary, this problem 
needs to be completely determined the data of the total mixture mass (say M m ) and of the total gas mass 
(say M g ) present in the computational domain. A natural way to impose these two conditions is to add to this 
problem two regularizing terms in the mass balance and the gas mass balance: 



c(h) \K\ 
c(h) \K\ 



g Zk) ~ W 



Mg 

ZK ~W\ 



E <k e p ' z (PK,ZK)-v- K Q v- z { PL ,z L ) = VKeM 

a=K\L 

z k ~ v- K z L = VK e M 



=K\L 



where c(h) is a regularization parameter tending to zero with the size of the mesh. In this case, the present 
existence theory directly applies, provided that the momentum balance equation remains linear with respect 
to the velocity. Of course, under the same restriction, this is true also for an implicit discretization of a 
time-dependent problem. 

In view of the stability results provided for the advection operator, adding such a term to the first relation 
of the problem (i.e. the momentum balance) should lead to a rather straightforward extension of the present 
existence result; the advection term would be multiplied by the homotopy parameter and the stability (i.e. an 
analogue to estimate (IbTj) ) would stem from the diffusion term. Note that, in this case, to keep the stability 
of the advection term, a regularization term consistent with the mass balance one should also be introduced in 
the momentum balance equation. 

We have shown in this paper that, with a Darcy's law for the drift velocity and a particular discretization for 
this term, the drift term is dissipative. So this term does not prevent to obtain stability estimates as (|6Tj) ; this 
suggests that the existence theory developped here may perhaps be extended to the complete drift flux model. 

Finally, we have not dealt in this study with the case where liquid monophasic zones (z = 0) exists in the flow. 
In such zones, the pressure changes of mathematical nature: it is no more a parameter entering the equation of 
state and determined by the local density, but a Lagrange multiplier for the incompressibility constraint. Note 
that this fact is already underlying in the present study: indeed, the incompressibility of the liquid prevents 
to derive L°° estimates for the pressure from L°° estimates for the density (which are readily obtained using 
a conservation argument), and we must invoke to this purpose the stability of the discrete gradient (i.e. the 
discrete inf-sup condition), that is typically the argument allowing to control the pressure in incompressible flow 
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problems. However, obtaining a priori estimates when z may vanish in the flow seems a difficult task, which 
should deserve more efforts. On the contrary, obtaining existence results for two barotropic phases seems to be 
rather simpler than the analysis performed here. 
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